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Abstract 

Isothermal  two-  and  three-dimensional  polymer  electrolyte  membrane  (PEM)  fuel  cell  cathode  flow  field  models  were  implemented  to  study 
the  behavior  of  reactant  and  reaction  product  gas  flow  in  a  parallel  channel  flow  field.  The  focus  was  on  the  flow  distribution  across  the  channels 
and  the  total  pressure  drop  across  the  flow  field.  The  effect  of  the  density  and  viscosity  variation  in  the  gas  resulting  from  the  composition  change 
due  to  cell  reactions  was  studied  and  the  models  were  solved  with  governing  equations  based  on  three  different  approximations.  The  focus  was  on 
showing  how  a  uniform  flow  profile  can  be  achieved  by  improving  an  existing  channel  design.  The  modeling  results  were  verified  experimentally. 
A  close  to  uniform  flow  distribution  was  achieved  in  a  parallel  channel  system. 

©  2006  Elsevier  B.Y.  All  rights  reserved. 
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1.  Introduction 

Fuel  cell  is  an  electrochemical  device  that  converts  the  chem¬ 
ical  energy  of  reactants  into  electricity  and  heat.  Fuel  cells 
typically  have  comparatively  high  efficiencies  and  energy  den¬ 
sities,  in  addition  to  which  they  have  potential  as  economically 
friendly  power  sources.  Their  properties  make  them  alternative 
power  sources  for  many  applications  ranging  from  portable  elec¬ 
tronics  and  vehicles  to  distributed  energy  production  and  power 
plants.  In  this  study  the  focus  is  on  the  polymer  electrolyte  mem¬ 
brane  fuel  cell  (PEMFC),  but  the  results  can  also  be  applied  to 
other  fuel  cell  types.  The  PEM  fuel  cell  operates  in  the  tem¬ 
perature  range  of  liquid  water,  though  higher  temperature  PEM 
fuel  cells  are  also  being  developed.  PEMFC  is  in  particular  suit¬ 
able  for  small-scale  applications  ranging  from  less  than  a  watt 
to  several  kilowatts. 

One  of  the  requirements  for  good  cell  performance  is  that  the 
reactants  must  be  distributed  as  uniformly  as  possible  across  the 
active  area  of  the  cell.  This  is  especially  important  on  the  cath¬ 
ode  side,  where  the  reaction  kinetics  is  comparatively  slow  and 
thus  forms  one  of  the  major  performance  limiting  factors  in  a 
PEMFC,  the  cathode  mass  transfer  overpotential.  A  non-uniform 


*  Corresponding  author.  Tel.:  +358  9  451  3209;  fax:  +358  9  451  3195. 
E-mail  address:  Suvi.Karvonen@hut.fi  (S.  Karvonen). 

0378-7753/$  -  see  front  matter  ©  2006  Elsevier  B.V.  All  rights  reserved, 
doi:  10. 1016/j  .jpowsour.2006.04. 145 


flow  distribution  results  in  a  non-uniform  reactant  distribution, 
leading  to  insufficient  amounts  of  reactants  on  some  areas  and 
inhomogeneous  current  production.  In  addition  to  having  an 
adverse  effect  on  cell  performance,  this  can  lead  to  tempera¬ 
ture  gradients  across  the  active  area  of  the  cell,  a  phenomenon 
that  in  an  extreme  case  may  damage  the  membrane.  The  flow  dis¬ 
tribution  properties  also  affect  the  water  removal  from  the  cell, 
and  thus  a  non-uniform  flow  distribution  can  cause  flooding  in 
the  cell. 

In  a  fuel  cell  stack,  the  reactant  flow  is  typically  directed  to 
each  unit  cell  with  a  component  known  as  the  flow  field  plate, 
which  also  functions  as  a  mechanical  support  structure  and  an 
electrical  contact.  The  flow  field  plate  directs  the  gas  flow  into  the 
gas  diffusion  layer  through  a  channel  system,  usually  molded, 
etched  or  machined  on  the  surface  of  the  plate.  The  flow  distribu¬ 
tion  in  the  channel  system  is  determined  by  the  relative  hydraulic 
resistances  of  the  channel  system.  The  most  common  channel 
configurations  are  the  parallel,  serpentine  and  interdigitated  con¬ 
figurations  and  their  combinations,  studied  in,  e.g.  [1-9].  The 
channel  system  can  also  be  replaced  with  a  porous  metal  net  or 
foam  plate,  see,  e.g.  [10].  The  parallel  channel  configuration, 
which  was  studied  in  this  work,  typically  has  a  small  hydraulic 
resistance  and  thus  does  not  generate  a  large  pressure  drop  across 
the  cell.  On  the  other  hand,  the  parallel  channel  flow  field  often 
has  a  non-uniform  flow  distribution  and  is  thus  more  susceptible 
to  flooding,  as  many  authors  have  concluded;  see,  e.g.  [5,6,9]. 
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Nomenclature 

A  area  (m2) 

D  diffusion  coefficient  (m2  s-1) 

Dh  hydraulic  diameter  (m) 

e  unit  charge  (1.6022  x  10“ 19  C) 

F  Faraday’s  constant  (96485  C  mol-1) 
g  gravitational  acceleration  (9.8067  m  s-2) 

i  current  density  (Am-2) 

L  characteristic  length  (m) 

M  molar  mass  (kg  mol- 1 ) 

n  surface  normal  vector 

N  number 

Na  Avogadro’s  constant  (6.022169  x  1023  mol-1) 
p  pressure  (Pa) 

po  atmospheric  pressure  (101.315  kPa) 

R  molar  gas  constant  (8.315  J  mol  1  K-1^ 

Re  Reynolds  number 

s  distance  (m) 

Scon  source  term  for  continuity  equation  (kg  m-3  s  - 1 ) 
Sns  source  term  for  Navier-Stokes  equation 

(kg  m  2  s-2) 

Sh  Sherwood  number 

t  surface  tangential  vector 

T  temperature  (K) 

u  total  velocity  of  the  fluid  (ms-1) 

u  fluid  velocity  vector  (ms-1) 

uq  inlet  velocity  (ms-1) 

U  characteristic  velocity  (ms-1) 

V  volume  (m3) 

Vm  molar  volume  in  343  K  (0.0278  m3  mol-1) 
v  molar  fraction 

z  number  of  electrons  involved  in  a  reaction 

Greek  symbols 

a  water  transport  coefficient 

A  surface  roughness  (m) 

£  porosity 

T]  dynamic  viscosity  (Pas) 

X  stoichiometric  constant,  2 

p  density  (kg  m- 3 ) 


Subscripts  and  superscripts 

act  active  area 

air  air 

atm  atmospheric 

ave  average 

Ar  argon 

eh  channel 

cath  cathode 

eff  effective 

H2O  water 

in  inlet 

lim  limiting 

max  maximum 


N2  nitrogen 

O2  oxygen 

react  reaction  participant 
sat  saturated  vapor 

tot  total 

v  vapor 


However,  these  problems  can  at  least  partially  be  avoided  with 
careful  design  of  the  flow  field  system.  A  uniform  flow  distri¬ 
bution  achieved  with  the  parallel  channel  system  is  presented  in 
this  work.  Consequently,  the  main  result  of  this  work  is  to  show 
by  example  that  uniform  flow  distributions  can  be  achieved  with 
parallel  channel  flow  fields  with  relatively  slight  changes  in  the 
flow  field  design. 

The  flow  distributions  and  pressure  losses  across  the  chan¬ 
nel  systems  were  studied  with  both  two-  and  three-dimensional 
one-phase  models  based  on  the  Navier-Stokes  and  continuity 
equations.  The  changes  in  the  gas  density  and  viscosity  along 
the  channels  that  result  from  the  cell  reactions  were  taken  into 
account  in  the  modeling.  The  cathode  distribution  was  studied 
according  to  three  different  approximations  and  the  correspond¬ 
ing  results  were  compared  in  order  to  find  out  the  error  induced 
by  each  approximation. 

The  modeling  data  and  experimental  results  showed  the  flow 
distribution  of  the  original  three-dimensional  parallel  channel 
system  to  be  polarized.  Based  on  the  results,  the  local  hydraulic 
resistances  of  the  channel  system  were  adjusted  through  modi¬ 
fication  of  the  gas  distributor  channel.  As  a  result,  the  modeled 
polarization  was  reduced  and  a  close  to  uniform  flow  profile 
was  achieved.  The  modeling  results  were  also  experimentally 
verified  and  the  experimental  results  were  in  agreement  with  the 
modeling  data. 


2.  Modeling 

2.7.  Navier-Stokes  and  continuity  equations 


The  modeling  domains  consisted  of  the  volume  in  3D  and 
cross-sectional  area  in  2D  of  the  modeled  channel  systems.  Other 
fuel  cell  components  such  as  the  gas  diffusion  layer  and  the 
membrane  electrode  assembly  were  excluded  from  the  model 
since  the  focus  in  this  work  was  on  the  flow  distribution  in  the 
flow  channels. 

The  incompressible  fluid  flow  is  governed  by  the  time- 
independent  Navier-Stokes  and  continuity  equations: 

-you  •  Vu  +  yOg  -  VjO  +  V  •  (^(Vu  +  VuT))  =  SnS  (1) 

V  .  (you)  =  Scon  (2) 


which  apply  to  a  laminar  flow.  The  flow  region,  laminar  or  tur¬ 
bulent,  is  defined  by  the  Reynolds  number,  Re: 


Re  = 


pu  D\\ 
0 


(3) 
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A  flow  is  laminar  when  Re  <  2000.  In  the  studied  flow  field  mod¬ 
els  the  Reynolds  number  is  below  160  and  Eqs.  (1)  and  (2)  apply 
to  the  problem  under  study. 

The  source  terms  in  Eqs.  (1)  and  (2),  Scon  and  Sns?  are  cal¬ 
culated  from  the  change  in  gas  composition  caused  by  the  cell 
reactions.  The  half  and  full  cell  reactions  of  the  PEMFC  are: 


anode  :  H2  — >  2H+  +  2e  (4) 

cathode  :  5O2  +  2H+  +  2e_  — >  H2O  (5) 

full  cell :  ±02  +  H2  -*  H20  (6) 


The  continuity  equation  source  term  S^on  is  calculated  from  the 
mass  difference  between  the  consumed  oxygen  and  produced 
water  molecules.  The  Navier-Stokes  source  term  Sns  is  cal¬ 
culated  based  on  the  assumption  that  the  kinetic  energy  and 
momentum  of  the  oxygen  consumed  in  the  reactions  are  lost 
into  the  gas  diffusion  layer  so  that  the  reaction  product  water 
has  no  initial  velocity.  Thus,  the  momentum  and  mass  source 
terms  can  be  written  as 


Sns 

5con 


iAact 

zFVch 


Mq2  u 


— — — — - (oi A/h2 o  -  M02) 
ZFV eh  2  2 


(7) 

(8) 


where  Vch  is  the  volume  of  the  channels  crossing  the  active  area 
of  the  cell.  However,  it  is  likely  that  some  of  the  kinetic  energy 
and  momentum  of  the  consumed  oxygen  is  in  fact  transferred 
to  the  produced  water  molecules,  but  estimating  the  magnitude 
of  this  phenomenon  is  difficult.  Therefore,  for  comparison,  the 
models  were  also  solved  using  an  alternative  Navier-Stokes 
source  term  calculated  by  assuming  that  the  kinetic  energy  and 
momentum  of  the  consumed  oxygen  molecules  are  transferred 
to  the  produced  water  molecules  so  that  their  initial  velocity 
equals  the  average  fluid  velocity  in  the  channel: 

Sns, 2  =  act  (aMu2o  Mo2)u  (9) 

ZP  Ech 


The  solutions  corresponding  to  the  different  Navier-Stokes 
source  terms  Sns  and  Sns, 2  gave  close  to  equal  results,  as 
the  relative  differences  in  the  flow  velocities  were  in  the  order 
of  10-3.  Thus,  the  possible  error  made  in  approximating  the 
Navier-Stokes  source  term  by  Eq.  (7)  should  be  negligible. 


2.2.  Approximations 


With  the  parameters  of  the  3D  model,  D h  =  0.67  mm  and 
Re max  =  160,  Eq.  (8)  gives  Z\iim  =  210  p,m.  This  is  clearly  larger 
than  the  surface  roughness  of  typical  flow  field  plate  materials 
such  as  graphite,  polymer  composite  and  steel,  whose  surface 
roughness  is  in  the  order  of  a  few  10  [xm  or  less,  see,  e.g.  [12,13]. 
However,  one  of  the  channel  walls  is  formed  by  the  gas  diffu¬ 
sion  layer  where  there  is  mass  transfer  through  the  surface  that 
should  be  taken  into  account.  Nevertheless,  including  this  effect 
would  complicate  the  modeling  and  consequently  was  excluded 
here. 

The  effect  of  gas  cross-over  between  channels  can  be 
neglected  since  the  Sherwood  number: 

Sh  =  >  where  Deff  =  Ds15  (12) 

is  in  the  order  of  1 03 ,  i.e.  the  flux  in  the  channels  is  three  orders  of 
magnitude  larger  than  the  diffusive  flux  in  the  gas  diffusion  layer. 
The  effect  of  possible  convective  flow  between  the  channels  is 
also  insignificant  since  the  pressure  differences  between  two 
parallel  channels  are  very  small,  at  largest  a  few  pascals,  which 
was  determined  from  the  solved  models. 

The  current  density  and  temperature  in  the  cell  are  assumed  to 
be  constant  across  the  active  area.  This  is  usually  not  the  case  in 
a  real  fuel  cell,  where  the  current  density,  temperature  and  flow 
distribution  are  all  interconnected  and  also  depend  on  external 
factors  such  as  the  cooling  system  of  the  cell.  Taking  all  these 
issues  into  account  would  have  made  the  models  very  complex 
and  as  a  consequence  required  a  lot  of  computing  capacity,  which 
was  the  reason  why  constant  values  for  these  parameters  were 
assumed  in  this  work. 

The  assumption  of  one-phase  flow,  i.e.  no  liquid  water  in  the 
channels,  is  justified  if  the  relative  humidity  remains  lower  than 
100%.  Taking  into  account  the  reaction  product  water  of  which 
half  is  assumed  to  leave  the  cell  through  the  cathode  side,  this 
applies  if  the  relative  humidity  of  the  inlet  gases  is  below  64%, 
which  was  calculated  assuming  a  stoichiometry  of  two  and  a  con¬ 
stant  cell  temperature  of  343  K.  These  values  correspond  to  the 
parameter  values  used  in  the  modeling.  In  many  real  fuel  cells, 
the  fluid  may  be  in  a  two-phase  flow,  and  the  modeling  results 
gained  here  do  not  necessarily  apply  in  these  cases.  However, 
high-temperature  and  low-humidity  membranes  that  function 
under  one-phase  flow  operating  conditions  such  as  assumed  in 
the  modeling  have  been  developed  (see,  e.g.  [14-17])  and  thus 
the  assumption  should  be  valid  for  several  real  fuel  cells. 


Eq.  (1)  can  be  simplified  for  modeling  purposes.  According  to 
dimensional  analysis  the  gravity  force  term  pg  can  be  excluded, 
since  its  weight  is  approximately  10-3  times  the  weight  of  the 
inertial  term: 

lpgl  *  n  =  *  10-3  (10) 

|pu  •  Vu|  pU2/L  U2 

The  channel  walls  in  a  flow  system  can  be  assumed  smooth,  if 
the  surface  roughness  of  the  wall  material  is  smaller  than  the 
limiting  surface  roughness  characteristic  to  that  system  [11]: 

Z\i,m  =  17.85  £>hfle-0'875  (1 1) 


2.3.  Model  properties 

A  schematic  of  the  2D  geometry  is  displayed  in  Fig.  1 .  The 
design  of  the  3D  geometry  is  similar  in  principle,  but  here  the 
distributor  channel  and  the  parallel  channels  crossing  the  cell 
are  in  different  planes  and  the  channels  are  divided  into  groups 
of  five  channels.  A  schematic  of  the  3D  geometry  is  presented 
in  Fig.  2.  The  height  of  the  distributor  channel  is  1  mm  and  the 
height  of  the  parallel  channels  is  0.5  mm.  The  cylinders  that 
connect  the  distributor  channel  and  the  parallel  channels  have  a 
radius  of  0.5  mm  and  height  of  2  mm.  Each  cylinder  distributes 
the  flow  to  five  parallel  channels  and  certain  periodicity  resulting 


S.  Karvonen  et  al.  /  Journal  of  Power  Sources  161  (2006)  876-884 


879 


Supporting  rib 


o 

OJ 


a> 

c 

c 

CD 

sz 

o 


3 

-Q 

00 

b 


Outlet 


Ir  let 


(7) 

< 

3 

3 

CD 


cr 

o 

c 

3 

Q. 

Q) 


126 


Fig.  1.  A  schematic  of  the  2D  geometry. 


from  this  can  also  be  seen  in  the  channel  flow  velocities  as  will 
be  discussed  later  in  Section  3.2.  At  the  inlet  boundaries  marked 
in  Figs.  1  and  2,  the  velocity  is  fixed:  u (w,  v,  w)  =  (wo,  0,  0).  At 
the  outlet  boundaries  also  marked  in  Figs.  1  and  2,  the  pressure 
is  fixed  to  zero,  p  =  0.  The  absolute  value  of  the  outlet  pressure 
does  not  affect  the  fluid  behavior  since  that  is  dependent  only 
on  the  pressure  gradients  within  the  flow  field,  which  can  be 
determined  from  Eq.  (1).  Due  to  symmetric  channel  geometries, 
only  one  half  of  the  cathode  flow  field  was  modeled.  The  sym¬ 
metry  boundaries  were  modeled  with  the  symmetry  boundary 


Outlet 


co 

*< 

3 

3 

CD 

CD 

Q_ 

CQ 

CD 


Fig.  2.  A  schematic  of  the  3D  geometry. 


condition: 


u  n  =  0  (13) 

ti  •  ?7(Vii  +  VuT)n  =  0  (14) 

and 

t2  •  ?7(Vu  +  VuT)n  =  0  (15) 


The  majority  of  the  boundaries,  corresponding  to  impermeable 
channel  walls,  were  governed  with  the  no-slip  condition  u  =  0. 
It  should  be  noted  that  the  no-slip  condition  was  applied  also 
to  the  wall  formed  by  the  gas  diffusion  layer  in  order  to  avoid 
further  complexity  in  the  modeling,  despite  the  fact  that  in  a  real 
fuel  cell  mass  transfer  exists  through  this  boundary. 

The  inlet  flow  velocity  is  calculated  from  the  current  density 
i  that  was  used  as  a  solver  parameter: 

«o?  RT  RT  AAact 
w0  =  — — - = - — - i  (16) 

Ao2,in  PO^in  -^02,in P0  ^cath^VA^in 

where  Aact  is  the  active  area  of  the  cell  and  Ain  the  cross-sectional 
area  of  the  inlet.  The  change  in  the  density  and  viscosity  of  the 
gas  in  the  channels  is  calculated  from  the  molar  fractions  of 
oxygen  and  water: 


,  x  ^02  ( wch^act  ^ac t  s 

X02(s)  =  —  =  — 7} - *02,in - “7 

^tot  V  £cathF  / 


X 


uc\\A 


Vr 


m 


act  iAact  s  ,  ~  *Aact  S 
- b  2a - 

^cathF  l  ^cathF  / 


iV- , 


—  (  x02,in 


m 


x  1  -  (2a  -  1) 


^ch^cathF  / 

iVm 


-1 


^ch^cathF  l 


(17) 


■^HoOW  = 


^h2o 

^tot 


=  *H20,in  +  2a 


iVi 


m 


x  1  -  (2a  -  1) 


iVr 


^ch^cathF  / 
-1 


m 


^ch^cathF  / 


(18) 


The  individual  channel  velocities  wc h  were  calculated  by  numer¬ 
ical  integration  across  the  channel  volume  separately  for  each 
channel: 


1  v- 

wch  —  —  /  UiVi  (19) 

Fch^ 

i= 1 

where  Vch  is  the  volume  (area  in  2D)  of  the  channel,  N  the 
number  of  calculation  points,  w;  and  Vi  are  the  velocity  at  point 
i  and  the  weight  factor  at  point  i  (volume  in  3D  and  area  in  2D 
of  the  space  represented  by  point  (/)),  respectively.  Thus,  wc h  is 
the  average  velocity  in  the  channel. 

In  the  modeling,  it  was  assumed  that  the  water  transport  coef¬ 
ficient  a  =  0.5,  i.e.  one  half  of  the  product  water  leaves  through 
the  anode  side.  Thus,  the  term  2a  —  1  =  0  and  the  equations  for 
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the  molar  fractions  simplify  to: 


*02  (^)  —  *02,  in 


iV- , 


m 


*H2oC?)  =  *H20,in  +  2a 


UchZcathF  l 

iVm 


UchZcathF  l 


(20) 

(21) 


These  equations  were  used  in  the  models  with  the  addition  that 
the  molar  fraction  of  oxygen  cannot  be  negative  and  the  molar 
fraction  of  water  has  a  maximum  value  corresponding  to  the 
situation  where  all  oxygen  has  been  consumed.  No  acceptable 
solution  can  exist  outside  these  limits. 

Using  the  ideal  gas  law  and  assuming  that  dry  air  is  a 
mixture  of  oxygen,  nitrogen  and  argon,  the  density  of  the  gas 
is  calculated  from: 


(Patm  +  P)(MH2oXu2o(s) 

+  Mo2Xo2(s )  +  Mn2IN2  +  MajXAt) 

RT 


(22) 


where  VN2and  xay  are  constant  when  a  =  0.5.  The  viscosity  of 
the  gas  mixture  is  [11]: 


^  +  *02  00  |  *n2  |  *Ar 


-l 


t?h2o 


V  o2 


0N2  OAr 


P(s) 

Pin 


(23) 


where  pm  is  the  density  of  the  dry  inlet  air.  The  density  and 
viscosity  values  of  the  relevant  gases  as  well  as  their  molar 
fractions  in  dry  standard  air  are  listed  in  Appendix  A. 

A  non-uniform  flow  profile  makes  it  easy  to  study  the  effect  of 
different  approximations  on  the  solution  since  the  differences  are 
usually  more  prominent  at  the  local  minima  and  maxima.  A  non- 
uniform  flow  profile  is  easier  to  accomplish  with  a  2D  geometry 
due  to  smaller  hydraulic  resistance,  and  thus  the  modeling  was 
done  both  in  two  and  three  dimensions.  In  both  dimensions,  three 
different  modeling  schemes  were  employed  to  study  the  effect 
of  the  density  and  viscosity  variation: 


1.  Constant  density  and  viscosity:  p  =  pair  and  rj  =  rjjn-  Zero 
source  terms.  Continuity  equation:  V-u  =  0. 

2.  Varying  density  and  viscosity:  p  =  p(s)  and  rj  =  rl(s)-  Zero 
source  terms.  Continuity  equation:  V  u  =  0. 

3.  Varying  density  and  viscosity:  p  =  p(s)  and  ri  =  rl(s)' 
Nonzero  source  terms:  5con  and  Sns-  Continuity  equation: 
V-(pu)  =  Scon. 


The  models  were  solved  with  the  commercially  available  par¬ 
tial  differential  equation  software  FEMLAB®.  The  calculations 
were  performed  over  a  64-bit  FEMLAB®  client-server  connec¬ 
tion.  The  server  computer  was  an  AMD  Athlon64  3500+  with 
4GB  RAM  and  40GB  of  swap-space.  The  operating  system  was 
SuSe  9.1  AMD64  Linux.  The  2D  geometry  was  modeled  with 
28  000  mesh  elements  resulting  in  160  000  degrees  of  freedom, 
whereas  the  3D  geometry  had  180  000  elements  and  1.2  million 
degrees  of  freedom.  In  each  element,  quadratic  Lagrange  poly¬ 
nomials  were  used  as  shape-functions  for  the  components  of 
the  velocity  field  while  linear  polynomials  were  used  for  pres¬ 
sure.  The  models  were  solved  to  as  high  current  densities  as 
the  FEMLAB®  solver  could  reach,  i.e.  0.35  A  cm-2  for  the  2D 


Fig.  3.  The  relative  channel  velocities  in  the  2D  model. 


geometry  and  0.4  A  cm-2  for  the  3D  geometry.  The  experimen¬ 
tal  parameters  corresponded  to  a  0.5  A  cm-2  current  density, 
which  is  the  designed  operating  current  for  the  studied  cell  and 
sufficiently  close  to  the  0.4  A  cm-2  current  density  of  the  3D 
model  for  the  flow  field  profiles  to  be  comparable. 

3.  Results 

3.1.  2D  geometry 

The  two-dimensional  modeling  domain  consisted  of  the 
cross-sectional  area  of  60  straight  parallel  channels  and  distrib¬ 
utor  channels  as  illustrated  in  Fig.  1.  The  results  discussed  here 
correspond  to  the  highest  current  density  at  which  the  model 
converged,  0.35  A  cm-2.  The  modeled  2D  flow  distribution  is 
illustrated  in  Fig.  3.  The  channel  velocities  are  taken  as  the  aver¬ 
age  velocities  in  the  channels  integrated  over  the  channel  region 
where  the  flow  is  stabilized  according  to  the  principle  that  was 
presented  in  Eq.  (19).  The  flow  distribution  is  strongly  polar¬ 
ized,  as  the  velocities  close  to  the  edges  of  the  flow  field  plate 
are  more  than  two  times  larger  than  the  smallest  channel  veloci¬ 
ties.  The  lines  corresponding  to  the  different  approximations  are 
close  to  indistinguishable.  Therefore,  the  differences  in  channel 
velocities  predicted  by  Schemes  1-3  are  illustrated  in  Fig.  4. 
Based  on  this  data,  the  maximum  differences  in  relative  individ- 


Fig.  4.  The  relative  differences  of  the  results  given  by  the  different  approxima¬ 
tions. 
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Max:  13.419 


Min:  -0.268 


Fig.  7.  The  relative  differences  between  the  different  modeling  schemes  of  the 
original  3D  geometry. 


Fig.  5.  The  pressure  distribution  in  the  2D  model.  The  pressure  dimension  is 
pascal  (Pa). 

ual  channel  velocities  between  Schemes  2  and  3  were  7.4%  and 
8.4%  for  Schemes  3  and  1.  The  total  pressure  losses  of  Schemes 
1  and  2  were  1.1%  and  2.3%  smaller  than  that  of  Scheme  3, 
respectively.  For  clarity,  the  pressure  distribution  across  the  2D 
geometry  corresponding  to  Scheme  3  with  0.35  A  cm-2  current 
density  is  illustrated  in  Fig.  5. 

3.2.  3D  geometry 

The  three-dimensional  model  corresponded  to  an  existing 
parallel  channel  flow  field  plate.  The  modeling  results  discussed 
here  are  the  solutions  corresponding  to  the  highest  current  den¬ 
sity  at  which  the  model  converged,  which  was  0.4  A  cm-2  for  the 
3D  geometry.  The  velocity  profiles  for  smaller  current  densities 
do  not  significantly  differ  from  the  ones  presented  here.  How¬ 
ever,  these  results  do  not  apply  to  significantly  higher  current 
densities  or  stoichiometric  ratios  where  the  increased  flow  rate 
causes  more  turbulence.  The  3D  flow-field  consisted  of  60  chan¬ 
nels  divided  into  12  five-channel  groups  and  distributor  channels 
such  as  illustrated  in  Fig.  2. 

The  relative  channel  velocities  are  presented  in  Fig.  6,  where 
a  distinctive  five-channel  periodicity  that  derives  from  the  chan¬ 
nel  system  design  can  be  seen  in  the  channel  velocities.  The 
flow  profile  of  the  3D  model  is  significantly  more  uniform  than 


Fig.  6.  The  relative  channel  flow  velocities  of  the  3D  model. 


that  of  the  2D  model  discussed  above.  The  largest  individual 
channel  velocity  is  16%  larger  than  the  smallest.  The  difference 
between  the  largest  and  smallest  five-channel  averages  is  12%. 
The  differences  in  relative  channel  velocities  between  Schemes 
1  and  3  are  illustrated  in  Fig.  7.  The  maximum  differences  are 
1.8%  (Scheme  3  versus  Scheme  1)  and  0.4%  (Scheme  3  versus 
Scheme  2),  significantly  smaller  than  those  of  the  2D  geometry. 
The  channel  velocities  in  Scheme  3  differ  more  from  the  veloci¬ 
ties  in  Schemes  1  and  2  in  the  channels  close  to  the  edges  of  the 
active  area.  This  follows  from  the  source  terms  in  Scheme  3  that 
take  into  account  the  momentum  that  is  lost  into  the  gas  diffusion 
layer.  The  total  gas  flow  sees  this  phenomenon  as  an  increase 
in  the  hydraulic  resistance  on  the  outlet  side.  Thus,  the  channel 
velocities  are  slightly  larger  on  the  outlet  side  in  Scheme  3.  In 
the  2D  model,  this  effect  is  more  difficult  to  see  since  the  total 
differences  between  Schemes  1  and  3  are  significantly  larger. 
The  total  pressure  losses  of  Schemes  1  and  2  were  9.0%  and 
7.3%  larger  than  that  of  Scheme  3,  respectively.  The  pressure 
distribution  across  the  modified  3D  geometry  corresponding  to 
Scheme  3  with  0.4  A  cm-2  current  density  is  illustrated  in  Fig.  8. 

One  of  the  objectives  of  this  work  was  to  achieve  a  parallel 
channel  flow  field,  where  the  flow  velocities  in  different  channels 


Max:  214 


Min:  0 


Fig.  8.  The  pressure  distribution  in  the  3D  model.  The  pressure  dimension  is 
pascal  (Pa)  and  the  length  dimension  is  meter  (m). 
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Fig.  9.  The  original  and  modified  distributor  channel  geometries. 


Fig.  10.  The  close  to  even  channel  velocity  profile  of  the  modified  3D  model. 

are  close  to  equal.  To  accomplish  this,  the  hydraulic  resistance  of 
the  inlet  distributor  channel  was  modified  by  changing  its  geo¬ 
metrical  shape.  The  modified  distributor  channel  was  optimized 
for  the  0-0.5  A  cm-2  current  density  with  stoichiometric  ratios 
between  1  and  2.  With  higher  flow  rates  the  velocity  distribution 
will  likely  become  more  polarized.  Schematics  of  the  original 
and  modified  inlet  distributor  channel  geometries  are  illustrated 
in  Fig.  9. 

The  velocity  profile  of  this  modified  3D  geometry  is  illus¬ 
trated  in  Fig.  10.  The  largest  individual  channel  velocity  is 
approximately  8%  larger  than  the  smallest  channel  velocity, 
which  means  that  the  1 6%  difference  of  the  original  3D  geometry 
has  been  reduced  to  half.  When  the  differences  in  each  five- 
channel  group  average  are  studied,  the  12%  value  of  the  original 
3D  geometry  has  diminished  to  3.7%.  Thus,  the  effect  of  the 
modification  is  significant.  The  differences  in  relative  channel 
velocities  between  Schemes  1  and  3  are  illustrated  in  Fig.  11. 
The  maximum  differences  are  2.1%  (Scheme  3  versus  Scheme 
1)  and  0.7%  (Scheme  3  versus  Scheme  2).  These  values  are 


Fig.  1 1 .  The  relative  differences  between  the  different  modeling  schemes  of  the 
modified  3D  geometry. 


Channel  number 

Fig.  12.  The  relative  velocities  of  the  original  3D  geometry  with  two  different 
meshes. 

slightly  larger  than  those  of  the  original  3D  geometry,  which 
implies  that  the  error  made  in  excluding  the  source  terms  or  the 
density  and  viscosity  variation  depends  on  the  model  geometry 
to  some  extent.  The  total  pressure  losses  of  Schemes  1  and  2 
were  8.7%  and  2.4%  larger  than  that  of  Scheme  3,  respectively. 

The  effect  of  the  mesh  size  on  the  results  of  the  3D  model  was 
studied  using  a  mesh  of  220  000  elements  for  the  3D  model  and 
comparing  the  results  to  the  modeling  results  presented  above 
with  180000  elements.  The  velocity  profiles  are  compared  in 
Fig.  12.  The  differences  are  small;  the  maximum  difference  in 
the  relative  channel  velocities  is  0.5%.  The  total  pressure  dif¬ 
ferences  with  the  different  meshes  are  within  0.2%.  The  effect 
of  the  mesh  was  also  studied  with  a  single  five-channel  group, 
where  the  number  of  elements  could  be  grown  to  four  times 
the  original.  The  relative  channel  velocities  were  still  within  the 
0.5%  error  marginal,  but  the  total  pressure  difference  grew  to 
5%.  This  gives  reason  to  expect  that  the  relative  channel  veloci¬ 
ties  should  be  fairly  reliable,  but  that  the  total  pressure  difference 
is  more  inaccurate.  Solving  the  3D  model  with  a  larger  number 
of  elements  was  not  realistic  since  with  the  available  capacity 
reaching  a  solution  for  the  model  took  from  15  to  40  h  with  the 
used  mesh. 

4.  Experimental  visualization 

The  modeling  data  of  the  original  and  modified  3D  geome¬ 
tries  were  verified  with  experiments.  The  experimental  proce¬ 
dure  was  simple:  dye  (water-soluble  black  ink)  was  mixed  to 
the  fluid  flow  in  pulses,  and  the  progress  of  the  dye  pulse  was 
recorded  with  a  digital  camera.  Since  mixing  dye  into  the  gas 
flow  would  have  been  complicated,  the  type  of  fluid  was  changed 
from  gaseous  (air)  to  liquid  (water).  In  order  to  conserve  the 
behavior  of  the  flow,  the  Reynolds  number  must  remain  con¬ 
stant  in  accordance  with  dimensional  analysis.  Following  from 
Eq.  (3)  the  product  of  density  and  velocity  divided  by  viscos¬ 
ity  must  remain  constant  since  the  hydraulic  diameter  is  not 
changed.  The  values  of  density  and  viscosity  for  air  and  water 
are  listed  in  Table  1 .  A  similar  visualization  study  utilizing  the 
laser-induced  fluorescence  method  has  previously  been  carried 
out  by  Barreras  et  al.  [9]. 
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Fig.  13.  The  experimental  arrangement  of  the  flow  visualization  study.  Because  of  symmetric  cell  structure,  the  dye  was  mixed  to  only  one  of  the  two  inlet  flows. 


The  Reynolds  number  for  an  air  flow  in  the  channels  with  a 
flow  velocity  M=l.lms_1  (corresponding  to  current  density  of 
i  =  0.5  A  cm-2  at  a  stoichiometry  of  two  for  which  the  studied 
cell  was  designed),  is  39.  The  flow  velocity  of  water  in  the  chan¬ 
nels  that  results  in  the  same  Reynolds  number  is  0.060ms_1, 
which  corresponds  to  an  inlet  flow  velocity  of  approximately 
200  ml  min-1. 

The  flow  field  plate  and  endplates  were  made  of  transpar¬ 
ent  polycarbonate.  The  usual  gasket  material  was  replaced  with 
PTFE.  The  water  flow  was  directed  into  the  flow  field  plate 
assembly  by  pressurizing  an  air-space  in  a  water  tank  with  a 
constant  200  ml  min-1  air  flow  to  provide  a  stabile  flow  in  the 
cell  assembly.  Circulating  water  through  the  system  caused  a 
part  of  the  channels  to  be  blocked  by  air  bubbles  due  to  the  high 
surface  tension  of  water.  Mixing  small  amounts  of  soap  into  the 
water  mostly  prevented  this  phenomenon.  The  dye  was  injected 
with  a  pipette  to  the  stabilized  water  flow  a  few  centimeters 
before  the  flow  entered  the  cell  assembly.  The  advancing  dye 
pulse  was  recorded  with  a  Sony®  DSC-F828  digital  camera. 
The  experimental  setup  is  illustrated  in  Fig.  13. 

The  progress  of  the  dye  pulse  in  each  five-channel  group  was 
studied  from  the  recorded  video  shots.  The  flow  velocities  in 
each  of  the  12  five-channel  groups  were  calculated  from  the 
time  spent  by  the  dye  pulse  to  move  through  the  length  of  the 
channels.  These  calculations  were  performed  with  eight  sets  of 
experimental  data  and  the  results  were  averaged  to  give  the  final 
values.  The  error  in  the  channel  velocities  was  taken  for  both 
geometries  as  the  maximum  average  deviation  of  the  individual 
channel  velocities  of  each  measurement  set,  4.3%  for  the  original 
geometry  and  5.6%  for  the  modified  geometry. 

The  average  velocities  of  the  12  five-channel  groups  are  com¬ 
pared  with  the  modeled  data  of  the  original  3D  geometry  in 
Fig.  14.  The  modeled  data,  to  which  the  comparisons  were  made, 

Table  1 

Channel  velocities  corresponding  to  constant  Reynolds  number  (Re  =  39)  and 
the  densities  and  viscosities  of  air  and  water 

p  (kg  m-3 )  rj  (Pas)  Inlet  volume  flow  (ml  min- 1 ) 

Air  1.031401  2.018-5  4104 

Water  1000  0.00103  217 


Fig.  14.  The  experimental  and  modeling  results  of  the  original  3D  geometry. 

is  that  of  Scheme  1  (constant  density  and  viscosity),  since  that 
is  closest  to  the  experimental  conditions,  where  the  density  and 
viscosity  of  the  water  remain  approximately  constant.  Within  the 
error  limits,  the  experimental  results  fit  the  model  data,  though 
the  experimental  velocity  profile  is  more  polarized.  The  stronger 
polarization  is  likely  due  to  some  systematic  error  affecting  the 
measurements  such  as  the  surface  tension  of  water.  Another  pos¬ 
sibility  is  that  imperfections  in  the  flow  field  geometry  could  also 
have  been  left  in  the  manufacturing  process. 

The  experiments  were  also  performed  with  the  modified  3D 
geometry.  The  results  are  illustrated  in  Fig.  15,  and  the  flow 
profile  is  more  uniform  than  in  Fig.  14  with  the  original  geometry 
as  predicted  by  the  modeling  data.  Thus,  the  experiments  confirm 


Fig.  15.  The  experimental  and  modeling  results  of  the  modified  3D  geometry. 
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that  the  modified  geometry  should  be  more  advantageous  for  use 
in  a  flow  field  plate  since  the  reactants  should  be  distributed  more 
uniformly  across  the  active  area  of  the  cell. 

5.  Summary  and  conclusions 

The  behavior  of  fluid  flow  in  an  isothermal  parallel  channel 
system  was  modeled  with  the  finite  element  method.  Three  mod¬ 
eling  schemes  based  on  different  approximations  were  employed 
in  the  modeling  and  the  results  achieved  by  these  schemes  were 
compared  to  each  other.  The  modeled  3D  parallel  channel  geom¬ 
etry  based  on  that  of  a  real  cell  was  modified  according  to 
modeling  data  and  the  improved  parallel  channel  system  had 
nearly  uniform  flow  distribution.  The  modeled  flow  profiles  were 
also  experimentally  verified.  The  experiments  were  carried  out 
by  recording  the  progress  of  a  dye  pulse  in  the  parallel  channels. 
The  distribution  of  the  measured  channel  velocities  was  in  good 
correlation  with  flow  distribution  predicted  by  the  modeling. 

It  was  discovered  that  neglecting  the  density  and  viscosity 
variation  caused  by  the  cell  reactions  caused  in  the  case  of  the  2D 
geometry  at  maximum  8%  differences  in  the  individual  channel 
velocities.  In  the  case  of  the  close  to  uniform  flow  profile  of  the 
3D  model  this  was  reduced  to  2%.  The  significance  of  this  error 
varied  between  the  used  geometries,  which  suggests  that  in  some 
cases  the  density  and  viscosity  variation  can  be  neglected,  but 
that  this  does  not  hold  generally. 

However,  the  results  imply  that  in  many  cases  the  effect  of 
excluding  the  cell  reactions  on  the  flow  profile  is  negligible  and 
thus  the  optimization  of  the  flow  field  channel  system  can  be 
done  separately  from  general  cell  optimization.  Uncoupling  the 
flow  field  channel  optimization  from  the  larger  cell  optimization 
problem  should  strongly  reduce  the  required  computing  capac¬ 
ity.  However,  the  real  non-isothermal  temperature  profile  of  the 
cell  depends  on  the  cooling  system  and  is  likely  to  have  some 
effect  on  the  flow  distribution,  offering  a  subject  for  further  stud¬ 
ies. 

Based  on  these  conclusions,  the  3D  parallel  channel  sys¬ 
tem  was  optimized  so  that  a  close  to  uniform  flow  profile  was 
achieved.  Thus,  it  has  been  demonstrated  in  this  work  that  one  of 
the  major  problems  in  using  the  parallel  channel  system  can  be 
overcome  with  careful  design  of  the  flow  field  plate.  This  makes 
the  parallel  channel  flow  field  a  promising  alternative  due  to  its 
typically  small  pressure  losses.  The  modeling  results  were  ver¬ 
ified  with  experiments  and  the  experimental  results  were  found 
to  be  in  agreement  with  the  modeling  data. 


Table  A.  1 


The  properties  of  dry  standard  air  used  in  the  modeling 


Component 

Molar  mass 
(g  mol- 1 ) 

Molar 

fraction  in 

standard 
air  (%) 

r]  at  343  K 
(Pas) 

p  at  343  K 
(kgm-3) 

Nitrogen 

28 

78 

1.97  x  10“5 

0.995 

Oxygen 

32 

21 

2.29  x  10“5 

1.137 

Argon 

40 

1 

2.60  x  10“5 

1.421 

Water 

18 

0 

1.15  x  10“5 

0.64 
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Appendix  A 

The  properties  of  dry  standard  air  is  shown  in  Table  A.l. 
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